#ifndef _MESH_
#define _MESH_

#include "Domain.h"
#include "Dofs.h"
#include "Grid.h"
#include <map>
#include <functional>
#include <Eigen/Sparse>
/**
 * @brief 定义域划分成网格Mesh,Mesh包含小网格Grid,每个Grid的类型由TemplateElement确定
 * 
 * @tparam int DIM 维度
 */
template <unsigned int DIM>
class Mesh
{
public:
    /**
     * @brief 默认构造函数
     */
    Mesh() = default;
    /**
     * @brief 含参构造函数,根据给定的定义域和划分段数来划分网格。
     * @param domain 定义域
     * @param seg 每个维度要划分的段数
     */
    Mesh(Domain<DIM> domain, vector<int> _seg);//只适用于结构网格
    /**
     * @brief 默认析构函数
     */
    virtual ~Mesh() = default;
    /**
     * @brief 给定模板单元的阶数，在每个Grid上划定自由度。
     * @details 仅实现二维一阶、二阶模板单元
     * 
     * @param _order_of_element 模板单元的阶数
     */
    virtual void distributeDofs(int _order_of_element);
    /**
     * @brief 返回自由度的总个数
     * @details _GridList.size()*getNumDofs();
     * @return 
     */
    virtual int getTotalNumDofs();
    /**
     * @brief 返回网格的总个数
     * @details _GridList.size()
     * @return [description]
     */
    virtual int getNumGrid();
    /**
     * @brief 返回存储Grid的vector
     * @details _GridList
     * @return [description]
     */
	virtual vector<Grid<DIM>>& getGridVector()
    {
        return gridList;
    }
    /**
     * @brief 返回序号为_index的Grid
     * 
     * @param _index Grid的序号
     * @return _GridList[i];
     */
	virtual Grid<DIM> &getGrid(int _index);
    /**
     * @brief   返回vector，存储边界上的自由度
     * @details 对所有grid,依据seg找出边界的grid,再确定对应边界的自由度下标是哪几个
     *
     * @param component 试图对边界分段，并标号
     * @return 存储边界上自由度的vector
     */
    virtual vector<Dofs<DIM>>& getBoundaryDofs(int component = 0);
    /**
     * @brief 返回map，存放边界上自由度的全局标号，和边界函数(仅Dirichlet)在对应自由度上的函数值。
     *
     * @param _boundary_function 边界函数
     * @param _boundary_dofs 边界上的自由度
     */
	virtual std::map<int,double> interpolateBoundaryValues(std::function<double(Dofs<DIM>&)>& _boundary_function, vector<Dofs<DIM>>& _boundary_dofs);
    /**
     * @brief 对刚度矩阵和右端项应用Dirichlet条件
     * 
     * @param system_matrix 刚度矩阵
     * @param system_rhs 右端项
     */
    virtual void applyBoundaryValues(std::map<int,double>& boundary_values,Eigen::SparseMatrix<double>& system_matrix,Eigen::VectorXd& system_rhs);
    //假定所有单元类型、阶数都相同


protected:
    vector<Grid<DIM>> gridList;/*< 存放所有的Grid*/
	vector<Dofs<DIM>> boundaryDofs;/*< 存放所有边界上的自由度*/
	vector<int> seg;/*< 存放各维度分段数*/

};

template <unsigned int DIM>
Mesh<DIM>::Mesh(Domain<DIM> domain, vector<int> _seg)
{
    assert(DIM==2 && _seg.size() == DIM);
	Point<2> p0 = domain.getVertex(0);
	Point<2> p2 = domain.getVertex(2);
	double xLen = (p2[0] - p0[0])/(1.0*_seg[0]);
	double yLen = (p2[1] - p0[1])/(1.0*_seg[1]);
	for(int j = 0; j < _seg[1]; ++j)
	{
		for(int i = 0; i < _seg[0]; ++i)
		{
			vector<Point<2>> pnts;
			pnts.push_back({p0[0]+i*xLen,p0[1]+j*yLen});
			pnts.push_back({p0[0]+(i+1)*xLen,p0[1]+j*yLen});
			pnts.push_back({p0[0]+(i+1)*xLen,p0[1]+(j+1)*yLen});
			pnts.push_back({p0[0]+i*xLen,p0[1]+(j+1)*yLen});
			Grid<2> grid(pnts);
			gridList.push_back(grid);
		}
	}
	seg = _seg;

};

template <unsigned int DIM>
int Mesh<DIM>::getNumGrid()
{
    return gridList.size();
};

template <unsigned int DIM>
Grid<DIM> &Mesh<DIM>::getGrid(int _index)
{
    return gridList[_index];
};

template <unsigned int DIM>
vector<Dofs<DIM>>&  Mesh<DIM>::getBoundaryDofs(int component)
{
    //对所有grid,依据seg找出边界的grid,再确定对应边界的自由度下标是哪几个
	static vector<Dofs<2>> v;
     for(int k = 0; k < seg[0]; k++)
          v.push_back(gridList[k].GetDof(0));
     for(int k = 0; k <= seg[0] * (seg[1] - 1); k = k + seg[0])
          v.push_back(gridList[k].GetDof(3));
     for(int k = seg[0] - 1; k < seg[0] * seg[1]; k = k + seg[0])
          v.push_back(gridList[k].GetDof(1));
     for(int k = seg[0] * (seg[1] - 1); k < seg[0] * seg[1]; k++)
       v.push_back(gridList[k].GetDof(2));
     v.resize(2 * (seg[0] + seg[1]));
     return v;

};
template <unsigned int DIM>
std::map<int,double> Mesh<DIM>::interpolateBoundaryValues(std::function<double (Dofs<DIM>&)> &_boundary_function, vector<Dofs<DIM> > &_boundary_dofs)
{
	std::map<int,double> boundary_values;
	for(Dofs<DIM>& dofs : _boundary_dofs)
	{
		int index = dofs.getGlobalIndex();
		double value = _boundary_function(dofs);
		boundary_values.insert({index,value});
	}
	return boundary_values;
}

template<unsigned int DIM>
void Mesh<DIM>::applyBoundaryValues(std::map<int,double>& boundary_values,Eigen::SparseMatrix<double>& system_matrix,Eigen::VectorXd& system_rhs)
{
	for(auto& it : boundary_values)
    {
        //边界自由度的全局编号。
		int k = it.first; 
        //对应自由的边界函数值。
		double bnd_value = it.second;
		system_rhs[k] = bnd_value * system_matrix.coeffRef(k,k);
		for(Eigen::SparseMatrix<double>::InnerIterator it(system_matrix,k);it;++it)
		{
			int row = it.row();
			if(row == k)
				continue;
			system_rhs[row] -= system_matrix.coeffRef(k,row)* bnd_value;
			system_matrix.coeffRef(k,row) = 0.0;
			system_matrix.coeffRef(row,k) = 0.0;
		}

    }
};
//假定所有单元类型、阶数都相同
template <unsigned int DIM>
void Mesh<DIM>::distributeDofs(int _order_of_element)
{
    //以下示范DIM=2，_order_of_element = 1 | 2;
    assert(DIM ==2);
    int global_idx = 0;
	int  k = 0;
    if(_order_of_element == 1)
    {
        vector<Dofs<DIM>> _dofs;
        for (Grid<DIM> &grid : gridList)
        {
            _dofs.clear();
            vector<Point<DIM>> p = grid.GetVertices();
            _dofs.push_back({{p[0][0], p[0][1]}, k + k / seg[0]});
            _dofs.push_back({{p[1][0], p[1][1]}, k + k / seg[0] + 1});
            _dofs.push_back({{p[2][0], p[2][1]}, k + k / seg[0] + seg[0] + 2});
            _dofs.push_back({{p[3][0], p[3][1]}, k + k /seg[0] + seg[0] + 1});
            grid.SetDofs(_dofs);
			k++;
        }

    }
    else if(_order_of_element == 2)
    {
        vector<Dofs<DIM>> _dofs;
        for (Grid<DIM> &grid : gridList)
        {
            _dofs.clear();
            vector<Point<DIM>> p = grid.GetVertices();
            _dofs.push_back({{p[0][0], p[0][1]}, global_idx++});
            _dofs.push_back({{(p[0][0] + p[1][0])/2.0, (p[0][1] + p[1][1])/2.0}, global_idx++});
            _dofs.push_back({{p[1][0], p[1][1]}, global_idx++});
            _dofs.push_back({{(p[1][0]+p[2][0])/2.0,(p[1][1] + p[2][1])/2.0}, global_idx++});
            _dofs.push_back({{p[2][0],p[2][1]}, global_idx++});
            _dofs.push_back({{(p[2][0]+p[3][0])/2.0,(p[2][1]+p[3][1]/2.0)}, global_idx++});
            _dofs.push_back({{p[3][0],p[3][1]}, global_idx++});
            _dofs.push_back({{(p[3][0]+p[0][0])/2.0,(p[3][1]+p[0][1])/2.0}, global_idx++});
            _dofs.push_back({{(p[0][0]+p[2][0])/2.0,(p[0][1]+p[2][1])/2.0}, global_idx++});
            grid.SetDofs(_dofs);
        }
    }
    else
    {
        std::cout <<"抱歉，干不了。\n"; 
    }

};

template <unsigned int DIM>
int Mesh<DIM>::getTotalNumDofs()
{
    int res = 1;
	for(int& n : seg)
	{
		res *=(n+1);//一次单元的自由度就是段数相乘
	}
    return res;
};


template<unsigned int DIM>
class UnstructedMesh: public Mesh<DIM>
{
public:
	virtual void readMesh();
	virtual void writeMesh();
	virtual void distributeDofs(int _order_of_element);

};

class EasyMesh :public  UnstructedMesh<2>
{

};

template <unsigned int DIM>
class StructedMesh : public Mesh<DIM>
{
public:
	void generateGrid(Domain<DIM> domain, vector<int> seg);
	virtual void distributeDofs(int _order_of_element);
};

#else
#endif
